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First all-sky upper limits from LIGO on the strength of periodic gravitational waves 

using the Hough transform 
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We perform a wide parameter space search for continuous gravitational waves over the whole sky 
and over a large range of values of the frequency and the first spin-down parameter. Our search 
method is based on the Hough transform, which is a semi-coherent, computationally efficient, and 
robust pattern recognition technique. We apply this technique to data from the second science run 
of the LIGO detectors and our final results are all-sky upper limits on the strength of gravitational 
waves emitted by unknown isolated spinning neutron stars on a set of narrow frequency bands in the 
range 200-400 Hz. The best upper limit on the gravitational wave strain amplitude that we obtain 
in this frequency range is 4.43 x 10 - 23 . 
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I. INTRODUCTION 

Continuous gravitational signals emitted by rotating 
neutron stars are promising sources for intcrferometric 
gravitational wave detectors such as GEO 600 Q, Q , 
the Laser Interferometer Gravitational Wave Observa- 
tory (LIGO) H0, TAMA 300 @ and VIRGO @. There 
are several physical mechanisms which might cause a 
neutron star to emit periodic gravitational waves. The 
main possibilities considered in the literature are (i) non- 
axisymmetric distortions of the solid part of the star 
011,11113, (ii) unstable r-modes in the fluid 0, [H Q2 , 
and (iii) free precession (or 'wobble') dill. The de- 
tectability of a signal depends on the detector sensitiv- 
ity, the intrinsic emission strength, the source distance 
and its orientation. If the source is not known, the de- 
tectability also depends on the available computational 
resources. For some search methods the detectability of 
a signal also depends on the source model used, but an 
all-sky wide-band search such as detailed here can detect 
any of the sources described above. 
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Previous searches for gravitational waves from rotat- 
ing neutron stars have been of two kinds. The first is 
a search targeting pulsars whose parameters are known 
through radio observations. These searches typically use 
matched filtering techniques and are not very computa- 
tionally expensive. Examples of such searches are [l5j 
and |l6j which targeted known radio pulsars, at twice 
the pulsar frequency, using data from the first and sec- 
ond science runs of the GEO 600 and LIGO detectors 
[l7j . No signals were detected and the end results were 
upper limits on the strength of the gravitational waves 
emitted by these pulsars and therefore on their cllipticity. 

The second kind of search looks for as yet undiscovered 
rotating neutron stars. An example of such a search is 
[Isj in which a two-day long data stretch from the Ex- 
plorer bar detector is used to perform an all sky search in 
a narrow frequency band around the resonant frequency 
of the detector. Another example is [l^ which uses data 
from the LIGO detectors to perform an all sky search in 
a wide frequency band using 10 hours of data. The same 
paper also describes a search for a gravitational wave sig- 
nal from the compact companion to Sco X-l in a large 
orbital parameter space using 6 hours of data. The key is- 
sue in these wide parameter space searches is that a fully 
coherent all-sky search over a large frequency band using 
a significant amount of data is computationally limited. 
This is because looking for weak continuous wave signals 
requires long observation times to build up sufficient sig- 
nal to noise ratio; the amplitude signal-to-noise ratio in- 
creases as the square root of the observation time. On the 
other hand, the number of templates that must be con- 
sidered, and therefore the computational requirements, 
scale much faster than linearly with the observation time. 
We therefore need methods which are sub-optimal but 
computationally less expensive [2£j, Wi, 122, 12J- I2H . Such 
methods typically involve semi-coherent combinations of 
the signal power in short stretches of data, and the Hough 
transform is an example of such a method. 

The Hough transform is a pattern recognition algo- 
rithm which was originally invented to analyze bubble 
chamber pictures from CERN [2^| . It was later patented 
by IBM |2jj, and it has found many applications in the 
analysis of digital images |2^. A detailed discussion of 
the Hough transform as applied to the search for con- 
tinuous gravitational waves can be found in [23l l24j . In 
this paper, we apply this technique to perform an all- 
sky search for isolated spinning neutron stars using two 
months of data collected in early 2003 from the second 
science run of the LIGO detectors (henceforth denoted 
as the S2 run). The main results of this paper are all-sky 
upper limits on a set of narrow frequency bands within 
the range 200-400 Hz and including one spin-down pa- 
rameter. 

Given the detector strain sensitivities during the S2 
run and its duration, it is unlikely to discover any neu- 
tron star, as described in the following Sec. II I II For this 
reason we focus here on setting upper limits. Substantial 
improvements in the detector noise have been achieved 
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since the S2 observations. A third science run (S3) took 
place at the end of 2003 and a fourth science run (S4) 
at the beginning of 2005. In these later runs LIGO in- 
struments collected data of improved sensitivity, but still 
less sensitive than the instruments' design goal. Several 
searches for various types of gravitational wave signals 
have been completed or are underway usin g d ata from 

the S2 and S3 runs EE E M M M El MM M we 

expect that the method presented here, applied as part 
of a hierarchical scheme and used on a much larger data 
set, will eventually enable the direct detection of periodic 
gravitational waves. 

This paper is organized as follows: Sec.[H]describes the 
second science run of the LIGO detectors; Sec. sum- 
marizes the current understanding of the astrophysical 
targets; Sec. IIVI reviews the waveform from an isolated 
spinning neutron star; Sec. presents the general idea 
of our search method, the Hough transform, and summa- 
rizes its statistical properties; Sec . I VII describes its imple- 
mentation and results on short Fourier transformed data. 
The upper limits are given in Sec. lVIII Sec. I Villi presents 
a validation of our search method using hardware injected 
signals, and finally Sec . IIXI concludes with a summary of 
our results and suggestions for further work. 



II. THE SECOND SCIENCE RUN 

The LIGO detector network consists of a 4 km inter- 
ferometer in Livingston, Louisiana (LI), and two inter- 
ferometers in Hanford, Washington, one 4 km and the 
other 2 km (HI and H2). Each detector is a power- 
recycled Michclson interferometer with long Fabry-Perot 
cavities in each of its orthogonal arms. These interfer- 
ometers are sensitive to quadrupolar oscillations in the 
space-time metric due to a passing gravitational wave, 
measuring directly the gravitational wave strain ampli- 
tude. 

The data analyzed in this paper were produced dur- 
ing LIGO's 59 day second science run. This run started 
on February 14 and ended April 14, 2003. Although the 
GEO detector was not operating at the time, all three 
LIGO detectors were functioning at a significantly bet- 
ter sensitivity than during SI, the first science run |l7j . 
and had displacement spectral amplitudes near 10 -18 m- 
Hz -1 / 2 between 200 Hz and 400 Hz. The strain sensitivi- 
ties in this science run were within an order of magnitude 
of the design sensitivity for the LIGO detectors. For a 
description of the detector configurations for S2 we refer 
the reader to Sec. IV and Sec. II. 

The reconstruction of the strain signal from the error 
signal of the feedback loop, used to control the differen- 
tial length of the interferometer arms, is referred to as 
the calibration. Changes in the calibration were tracked 
by injecting continuous, fixed-amplitude sinusoidal exci- 
tations into the end test mass control systems, and moni- 
toring the amplitude of these signals at the measurement 
error point. Calibration uncertainties at the three LIGO 




Frequency (Hz) 

FIG. 1: Characteristic amplitude detectable from a known 
generic source with a 1% false alarm rate and 10% false dis- 
missal rate, as given by Eq. 1171 . All curves use typical sensi- 
tivities of the three LIGO detectors during S2 and observation 
times corresponding to the up-time of the detectors during S2. 
The thin line is the expected characteristic amplitude for the 
same false alarm and false dismissal rates, but using the initial 
LIGO design goal for the 4 km instruments and an effective 
observation time of 1 yr. 



detectors during S2 were estimated to be smaller than 

11% m. 

The data were acquired and digitized at a rate of 
16384 Hz. The duty cycle for the interferometers, de- 
fined as the fraction of the total run time when the in- 
terferometer was locked (i.e., all interferometer control 
servos operating in their linear regime) and in its low 
noise configuration, were similar to those of the previous 
science run, approximately 37% for LI, 74% for HI and 
58% for H2. The longest continuous locked stretch for 
any interferometer during S2 was 66 hours for HI. The 
smaller duty cycle for LI was due to anthropogenic diur- 
nal low-frequency seismic noise which prevented opera- 
tions during the day on weekdays. Recently installed ac- 
tive feedback seismic isolation has successfully addressed 
this problem. 

Fig. n shows the expected sensitivity for the Hough 
search by the three LIGO detectors during S2. Those ho 
values correspond to the amplitudes detectable from a 
generic continuous gravitational wave source, if we were 
performing a targeted search, with a 1% false alarm rate 
and 10% false dismissal rate, as given by Eq. (|17fl . The 
differences among the three interferometers reflect dif- 
ferences in the operating parameters, hardware imple- 
mentation of the three instruments, and the duty cycles. 
Fig. n also shows the expected sensitivity (at the same 
false alarm and false dismissal rates) for initial LIGO 
4 km interferometers running at design sensitivity assum- 
ing an observation time of 1 year. These false alarm 
and false dismissal values arc chosen in agreement with 
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|l5| and |19| only for comparison purposes. Because of 
the large parameter search we perform here, it would be 
more meaningful to consider a lower false alarm rate, say 
of 10~ 10 and then the sensitivity for a targeted search 
would get worse by a factor 1.5. The search described in 
this paper is not targeted and this degrades the sensitiv- 
ity even further. This will be discussed in Sec. IVIII 

At the end of the S2 run, two fake artificial pulsar sig- 
nals were injected for a 12 hour period into all three LIGO 
interferometers. These hardware injections were done by 
modulating the mirror positions via the actuation control 
signals. These injections were designed to give an end- 
to-end validations of the search pipelines starting from 
as far up the observing chain as possible. See Sec. IVIIII 
for details. 



III. ASTROPHYSICAL TARGETS 

The target population of this search consists of isolated 
rotating neutron stars that are not observed in electro- 
magnetic waves. Current models of stellar evolution sug- 
gest that our Galaxy contains of order 10 9 neutron stars 
and that of order 10 5 are active pulsars [33]. Up to now, 
only of order 10 3 objects have been identified as neu- 
tron stars, either by observation as pulsars, or through 
their x-ray emission, of which about 90% are isolated 
HE HE HE [13 ■ Most neutron stars will remain unob- 
served electromagnetically for many reasons such as the 
non-pulsed emission being faint or the pulses being emit- 
ted in a beam which does not sweep across the Earth. 
Therefore, there are many more neutron stars in the tar- 
get population than have already been observed. 

Although there is great uncertainty in the physics of 
the emission mechanism and the strength of an individ- 
ual source, we can argue for a robust upper limit on the 
strength of the strongest source in the galactic population 
that is almost independent of individual source physics. 
The details of the argument and an overview of emission 
mechanisms can be found in a forthcoming paper |l9j |. 
Here we do not repeat the details but merely summa- 
rize the result. For an upper limit we make optimistic 
assumptions — that neutron stars are born rapidly rotat- 
ing and spinning down due to gravitational waves, and 
that they are distributed uniformly throughout the galac- 
tic disc — and the plausible assumption that the overall 
galactic birthrate 1 /tj, is steady. By converting these as- 
sumptions to a distribution of neutron stars with respect 
to gravitational wave strain and frequency, we find there 
is a 50% chance that the strongest signal between fre- 
quencies / min and / max has an amplitude of at least 



h « 4 x 10~ 



30yr\ /, 



n 



In- 



1/2 



(1) 



Of course, with less optimistic assumptions this value 
would be smaller. 

Comparing Eq. to Fig. ^ a search of S2 data is 
not expected to result in a discovery. However, it is 



still possible that the closest neutron star is closer than 
the typical distance expected from a random distribution 
of supernovae (for example due to recent star formation 
in the Gould belt as considered in Ref. 01 )■ ^ is also 
possible that a "blind" search of this sort may discover 
some previously unknown class of compact objects not 
born in supernovae. More importantly, future searches 
for previously undiscovered rotating neutron stars using 
the methods presented here will be much more sensitive. 
The goal of initial LIGO is to take a year of data at de- 
sign sensitivity, which means a factor 10 decrease in the 
amplitude strain noise relative to S2, and a factor 10 in- 
crease in the length of the data set. These combine to 
reduce ho to somewhat below the value in Eq. and 
thus initial LIGO at full sensitivity will have some chance 
of observing a periodic signal. 



IV. THE EXPECTED WAVEFORM 

In order to describe the expected signal waveform we 
will use the same notation as [l5j . We will briefly summa- 
rize it in the next paragraphs for convenience. The form 
of the gravitational wave emitted by an isolated spinning 
neutron star, as seen by a gravitational wave detector, is 



h(t) = F+(t, ^)h+{t) + F x (t, r/,)h x (t) 



(2) 



where t is time in the detector frame, i\) is the polariza- 
tion angle of the wave, and F + _ x are the detector an- 
tenna pattern functions for the two polarizations. If we 
assume the emission mechanism is due to deviations of 
the pulsar's shape from perfect axial symmetry, then the 
gravitational waves are emitted at a frequency which is 
exactly twice the rotational rate f r . Under this assump- 
tion, the waveforms for the two polarizations x are 
given by: 



h - 



1 



cos$(i) . 



h x = ho cost sin <!>(i) , 



(3) 
(4) 



where i is the angle between the neutron star's spin axis 
and the direction of propagation of the waves, and ho is 
the amplitude: 



h 



Wn 2 GI zz ef? 



(5) 



where G is Newton's gravitational constant, c the speed 
of light, I zz is the principal moment with the z-axis be- 
ing its spin axis, e := (I xx — I yy )/I zz is the equatorial 
cllipticity of the star, and d is the distance to the star. 

The phase takes its simplest form in the Solar Sys- 
tem Barycenter (SSB) frame where it can be expanded 
in a Taylor series up to second order: 

*(t) = $ + 2tt (f (T - To) + \j{T - Tof \ . (6) 



5 



Here T is time in the SSB frame and To is a fiducial start 
time. The phase $01 frequency /o and spin-down parame- 
ter / arc defined at this fiducial start time. In this paper, 
we include only one spin-down parameter in our search; 
as we shall see later in Sec. IVI Bl our frequency resolu- 
tion is too coarse for the higher spin-down parameters to 
have any significant effect on the frequency evolution of 
the signal (for the spin-down ages we consider). 

Modulo relativistic effects which are unimportant for 
this search, the relation between the time of arrival T of 
the wave in the SSB frame and in the detector frame t is 



where n is the unit vector from the detector to the neu- 
tron star, and r is the detector position in the SSB frame. 

The instantaneous frequency f{t) of the wave as ob- 
served by the detector is given, to a very good approxi- 
mation, by the familiar non-rclativistic Dopplcr formula: 

/(*)-/(*) = /(*) 1 ^- E , (8) 
c 

where f(t) is the instantaneous signal frequency in the 
SSB frame at time t: 

/(t) = /o + /(t-fc + ^^), 0) 

where to is the fiducial detector time at the start of the 
observation and Ar(t) = r(t) — r(to). It is easy to see 
that the Ar • n/c term can safely be ignored so that, to 
an excellent approximation 

f(t) = fo + f(t-to). (10) 

V. THE HOUGH TRANSFORM 

In this paper, we use the Hough transform to find 
the pattern produced by the Doppler shift ijSJl and the 
spin-down IjlUf) of a gravitational wave signal in the time- 
frequency plane of our data. This pattern is independent 
of the source model used and therefore of the emission 
mechanisms. We only assume that the gravitational wave 
signal is emitted by an isolated spinning neutron star. 

The starting point for our search is a set of data seg- 
ments, each corresponding to a time interval T coh . Each 
of these data segments is Fourier transformed to produce 
a set of N short time-baseline Fourier transforms (SFTs). 
From this set of SFTs, calculating the periodograms (the 
square modulus of the Fourier transform) and selecting 
frequency bins (peaks) above a certain threshold, we ob- 
tain a time-frequency map of our data. In the absence 
of a signal the peaks in the time-frequency plane arc dis- 
tributed in a random way; if signal is present, with high 
enough signal to noise ratio, some of these peaks will be 



distributed along the trajectory of the received frequency 
of the signal. 

The Hough transform maps points of the time- 
frequency plane into the space of the source parameters 
(/o,/,n). The result of the Hough transform is a his- 
togram, i.e., a collection of integer numbers, each repre- 
senting the detection statistic for each point in parame- 
ter space. We shall refer to these integers as the number 
count. The number counts are computed in the following 
way: For each selected bin in the SFTs, we find which 
points in parameter space are consistent with it, accord- 
ing to Eq. JSJ), and the number count in all such points is 
increased by unity. This is repeated for all the selected 
bins in all the SFTs to obtain the final histogram. 

To illustrate this, let us assume the source parameters 
are only the coordinates of the source in the sky, and this 
source is emitting a signal at a frequency f$. Moreover 
we assume that at a given time t a peak at frequency / 
has been selected in the corresponding SFT. The Hough 
transform maps this peak into the loci of points, on the 
celestial sphere, where a source emitting a signal with 
frequency fo could be located in order in order to produce 
at the detector a peak at /. By repeating this for all 
the selected peaks in our data we will obtain the final 
Hough map. If the peaks in the time-frequency plane 
were due only to signal, all the corresponding loci would 
intersect in a region of the Hough map identifying the 
source position. 

An advantage of the Hough transform is that a large 
region in parameter space can be analyzed in a single 
pass. By dropping the amplitude information of the se- 
lected peaks, the Hough search is expected to be com- 
putationally efficient, but at the cost of being somewhat 
less sensitive than others semi-coherent methods, e.g., 
the stack-slide search [2l|. On the other hand, discard- 
ing this extra information makes the Hough transform 
more robust against transient spectral disturbances be- 
cause no matter how large a spectral disturbance is in 
a single SFT, it will contribute at the most +1 to the 
number count. This is not surprising since the optimal 
statistic for the detection of weak signals in the presence 
of a Gaussian background with large non Gaussian out- 
liers is effectively cut off above some value 01 • This 
is, in practice, what the Hough transform does to large 
spectral outliers. 

With the above short summary at hand, we now give 
the relevant notation and equations that will be used 
later. For further details and derivations of the equations 
below, we refer the reader to p^ |. 

Frequency bins are selected by setting a threshold p th 
on the normalized power pk defined as 

2|x fc | 2 

Pk - 77, q 77 \ ) I 11 ) 

where Xk is the discrete Fourier transform of the data, the 
frequency index fc corresponds to a physical frequency of 
fk = ^'/Tohj and S n (fk) is the single sided power spectral 
density of the detector noise. The k th frequency bin is 
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selected if pk > p t h, and rejected otherwise. In this way, 
each SFT is replaced by a collection of zeros and ones 
called a peak-gram. 

Let n be the number count at a point in parameter 
space, obtained after performing the Hough transform 
on our data. Let p(n) be the probability distribution of 
n in the absence of a signal, and p(n\h) the distribution in 
the presence of a signal h(t). It is clear that < n < N, 
and it can be shown that for stationary Gaussian noise, 
p{n) is a binomial distribution with mean Nq where q is 
the probability that any frequency bin is selected: 

p(n)=(^q n (l-q) N - n . (12) 

For Gaussian noise in the absence of a signal, it is easy to 
show that pk follows an exponential distribution so that 
q = e~ Pth . In the presence of a signal, the distribution 
is ideally also a binomial but with a slightly larger mean 
Nr) where, for weak signals, ry is given by 

11 = q [l + P -^X + 0{\ 2 )} . (13) 

A is the signal to noise ratio within a single SFT, and for 
the case when there is no mismatch between the signal 
and the template: 

with h(f) being the Fourier transform of the signal h(t). 
The approximation that the distribution in the presence 
of a signal is binomial breaks down for reasonably strong 
signals. This is due to possible non-stationarities in the 
noise, and the amplitude modulation of the signal which 
causes 77 to vary from one SFT to another. 

Candidates in parameter space are selected by setting 
a threshold rt th on the number count. The false alarm 
and false dismissal rates for this threshold are defined 
respectively in the usual way: 

N nth-1 

a = P(n), P= Yl P{n\h). (15) 

n— n th n—Q 

We choose the thresholds (n th , p th ) based on the Neyman- 
Pearson criterion of minimizing f3 for a given value of a. 
It can be shown [2j| that this criteria leads, in the case 
of weak signals, large N , and Gaussian stationary noise, 
to p th w 1.6. This corresponds to q f=a 0.20, i.e., we select 
about 20% of the frequency bins from each SFT. This 
value of p th turns out to be independent of the choice of a 
and signal strength. Furthermore, rt th is also independent 
of the signal strength and is given by: 

n tb = Nq+ yj2Nq{l-q) erfc" 1 (2a) , (16) 

where erfc -1 is the inverse of the complementary error 
function. These values of the thresholds lead to a false 



dismissal rate (3 which is given in 23]. The value of 
(3 of course depends on the signal strength, and on the 
average, the weakest signal which will cross the above 
thresholds at a false alarm rate a and false dismissal (3 
is given by 

h ° = 5M ^^ (17) 

where 

5 = crfc _1 (2a)+crfc" 1 (2/3). (18) 

Equation (|17f) gives the smallest signal which can be de- 
tected by the search, and is therefore a measure of the 
sensitivity of the search. 



VI. THE SEARCH 
A. The SFT data 

The input data to our search is a collection of cal- 
ibrated SFTs with a time baseline T coh of 30 minutes. 
While a larger value of T coh leads to better sensitivity, 
this time baseline can not be made arbitrarily large be- 
cause of the frequency drift caused by the Dopplcr effect 
(and also the spin-down); we would like the signal power 
of a putative signal to be concentrated in less than half 
the frequency resolution l/T coh . It is shown in j2^ that 
at 300 Hz, we could ideally choose T coh up to ~ 60min. 
On the other hand, we should be able to find a significant 
number of such data stretches during which the interfer- 
ometers are in lock, the noise is stationary, and the data 
are labeled satisfactory according to certain data quality 
requirements. Given the duty cycles of the interferome- 
ters during S2 and the non-stationarity of the noise floor, 
it turns out that T coh = 30 min is a good compromise 
which satisfies these constraints. By demanding the data 
in each 30 minutes stretch to be continuous (although 
there could be gaps in between the SFTs) the number N 
of SFTs available for LI data is 687, 1761 for HI and 1384 
for H2, reducing the nominal duty cycle for this search. 

The SFT data are calibrated in the frequency domain 
by constructing a response function R(f, t) that acts on 
the error signal of the feedback loop used to control the 
differential length of the interferometer arms. The re- 
sponse function R(f, t) varies in time, primarily due to 
changes in the amount of light stored in the Fabry-Perot 
cavities of the interferometers. During S2, changes in 
the response were computed every 60 seconds, and an 
averaging procedure was used to estimate the response 
function used on each SFT. The SFTs are windowed and 
high-pass filtered as described in Sec. IV C 1 of 0]. No 
further data conditioning is applied, although the data 
are known to contain many spectral disturbances, includ- 
ing the 60 Hz power line harmonics and the thermally 
excited violin modes of test mass suspension wires. 
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B. The parameter space 

This section describes the portion of parameter space 
(/o, /,n) we search over, and the resolution of our grid. 
Our template grid is not based on a metric calculation 
(as in e.g. |2fl|21|), but rather on a cubic grid which cov- 
ers the parameter space as described below. Particular 
features of this grid are used to increase computational 
efficiency as described in Sec. IVI CI 

We analyze the full data set from the S2 run with a 
total observation time T obs ~ 5.f x 10 6 sec. The exact 
value of T obs is different for the three LIGO interferom- 
eters ■ We search for isolated neutron star signals in 
the frequency range 200-400 Hz with a frequency resolu- 
tion 



Sf=^— = 5.556 x 



10 _4 Hz. 



(19) 



The choice of the range 200-400 Hz for the analysis is 
motivated by the low noise level, and therefore our ability 
to set the best upper limits for ho, as seen from Fig. Q 

The resolution Sf in the space of first spin-down pa- 
rameters is given by the smallest value of / for which the 
intrinsic signal frequency docs not drift by more than a 
single frequency bin during the total observation time 



Sf = 



Sf_ 

T ohs 



1 



^obs^coh 



1.1 x 10- 10 Hz-s _1 



(20) 



We choose the range of values — / max < / < 0, where 
the largest spin-down parameter / max is about 1.1 x 
10 -9 Hz-s~ . This yields eleven spin-down values for each 
intrinsic frequency. In other words, we look for neutron 
stars whose spin-down age is at least r min = /// max . This 
corresponds to a minimum spin-down age of 5.75 x 10 3 yr 
at 200 Hz, and 1.15 x 10 4 yr at 400 Hz. These values 
of / max and r min are such that all known pulsars have a 
smaller spin-down rate than / max and, except for a few su- 
pernova remnants, all of them have a spin-down age sig- 
nificantly greater than the numbers quoted above. With 
these values of r min , it is easy to see that the second spin- 
down parameter can be safely neglected; it would take 
about 10 yr for the largest second spindown parameter 
to cause a frequency drift of half a frequency bin. 

As described in j23|, for every given time, value of the 
intrinsic frequency fo and spin-down /, the set of sky- 
locations n consistent with a selected frequency f(t) cor- 
responds to a constant value of v • n given by © ■ This is 
a circle in the celestial sphere. It can be shown that every 
frequency bin of width Sf corresponds to an annulus on 
the celestial sphere whose width is at least 



(66)» 



c6J_ 
v f 



(21) 



with v being the magnitude of the average velocity of the 
detector in the SSB frame. 



The resolution 50 in sky positions is chosen to be fre- 
quency dependent, being at most 58 = ^(58) mia . To 
choose the template spacing only, we use a constant value 
of v/c equal to 1.06 x 10~ 4 . This yields: 



9.3 x 10~ 3 rad x 



300 Hz 
/ 



(22) 



This resolution corresponds to approximately 1.5 x 10 5 
sky locations for the whole sky at 300 Hz. For that, we 
break up the sky into 23 sky-patches of roughly equal 
area and, by means of the stereographic projection, we 
map each portion to a plane, and set a uniform grid with 
spacing 50 in this stereographic plane. The stereographic 
projection maps circles in the celestial sphere to circles 
in the plane thereby mapping the annuli in the celestial 
sphere, described earlier, to annuli in the stereographic 
plane. We ensure that the dimensions of each sky-patch 
are sufficiently small so that the distortions produced by 
the stereographic projection are not significant. This is 
important to ensure that the number of points needed to 
cover the the full sky is not much larger than if we were 
using exactly the frequency resolution given by Eq. Q22p. 

This adds up to a total number of templates per 1 Hz 
band at 200 Hz ~ 1.9 x 10 9 while it increases up to 7.5 x 
10 9 at 400 Hz. 



C. The implementation of the Hough transform 

This section describes in more detail the implemen- 
tation of the search pipeline which was summarized in 
Sec. [V] The first step in this semi-coherent Hough search 
is to select frequency bins from the SFTs and construct 
the peak-grams. As mentioned in Sec.[V] our criteria for 
selecting frequency bins is to set a threshold of 1.6 on 
the normalized power l|ll|) . thereby selecting about 20% 
of the frequency bins in every SFT. 

The power spectral density S n appearing in Eq. (jll(l is 
estimated by means of a running median applied to the 
periodogram of each individual SFT. The window size we 
employ for the running median is w = 101 correspond- 
ing to 0.056 Hz [57^ . The runnin g m edian is a robust 
method to estimate the noise floor [43, |43, |4£] which has 
the virtue of discarding outliers which appear in a small 
number of bins, thereby providing an accurate estimate 
of the noise floor in the presence of spectral disturbances 
and possible signals. The use of the median (instead of 
the mean) to estimate the power spectral density intro- 
duces a minor technical complication (see appendixlAlfor 
further details). 

The next step is to choose a tiling of the sky. As 
described before, we break up the sky into 23 patches, 
of roughly equal area. By means of the stereographic 
projection, we map each portion to a two dimension 
plane and set a uniform grid with a resolution 58 in this 
plane. All of our calculations arc performed on this stere- 
ographic plane, and are finally projected back on to the 
celestial sphere. 
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FIG. 2: A partial Hough map (PHM) is a histogram in the 
space of sky locations obtained by performing the Hough 
transform on a single SFT and for a given value of the instan- 
taneous frequency. A total Hough map is obtained by sum- 
ming over the appropriate PHMs. The PHMs to be summed 
over are determined by the choice of spin-down parameters 
which give a trajectory in the time-frequency plane. 



In our implementation of the Hough transform, we 
treat sky-positions separately from frequencies and spin- 
downs. In particular, we do not obtain the Hough his- 
togram over the entire parameter space in one go, but 
rather for a given sky-patch, a search frequency fo and a 
spin-down / value. These are the so-called Hough maps 
(HMs). Repeating this for every set of frequency and 
spin-down parameters and the different sky-patches we 
wish to search over, we obtain a number of HMs. The 
collection of all these HMs represent our final histogram 
in parameter space. 

The HMs could be produced by using a "brute force" 
approach, i.e., using all the peaks in the time- frequency 
plane. But there is an alternative way of constructing 
them. Let us define a partial Hough map (PHM) as be- 
ing a Hough histogram, in the space of sky locations, 
obtained by performing the Hough transform using the 
peaks from a single SFT and for a single value of the 
intrinsic signal frequency and no spin-down. This PHM 
therefore consists of only zeros and ones, i.e., the col- 
lection of the annuli corresponding to all peaks present 
in a single peak-gram. Then each HM can be obtained 
by summing the appropriate PHMs produced from dif- 
ferent SFTs. If we add PHMs constructed by using the 
same intrinsic frequency, then the resulting HM refers to 
the same intrinsic frequency and no spin-down. But note 
that the effect of a spin-down in the signal is the same 
as having a time varying intrinsic frequency. This sug- 
gests a strategy to re-use PHMs computed for different 
frequencies at different times in order to compute the HM 
for a non-zero spin-down case. 

Given the set of PHMs, the HM for a given search fre- 
quency fo and a given spin-down / is obtained as follows: 
using Eq. I|10|) calculate the trajectory f(t) in the time- 
frequency plane corresponding to fo and /. If the mid 
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FIG. 3: The schematic of the analysis pipeline. The input 
data are the SFTs and the search parameters. The first step 
is to select frequency bins from the SFTs and generate the 
peak-grams. Then, the Hough transform is computed for the 
different sky patches, frequencies and spin-down values, thus 
producing the different Hough maps. The search uses LUTs 
that are computed for a given tiling of the sky-patch. The 
sky grid is frequency dependent, but it is fixed for the fre- 
quency range in which the LUTs are valid. Then, a collec- 
tion of PHMDs is built, and for each search frequency fo and 
given spin-down /, the trajectory in the time-frequency plane 
is computed and the Hough map obtained by summing and 
integrating the corresponding PHMDs. The code loops over 
frequency and spin-down parameters, updating the sky grid 
and LUTs whenever required. Statistical analyses are per- 
formed on each map in order to reduce the output size. 



time stamps of the SFTs are {ti} (i — 1 . . . N), calculate 
f(U) and find the frequency bin that it lies in; select the 
PHM corresponding to this frequency bin. Finally, add 
all the selected PHMs to obtain the Hough map. This 
procedure is shown in Fig. [21 

This approach saves computations because it recog- 
nizes that the same sky locations can refer to different 
values of frequency and spin-down, and avoids having to 
re-determine such sky locations more than once. Another 
advantage of proceeding in this fashion is that we can use 
look up tables (LUTs) to construct the PHMs. The basic 
problem to construct the PHMs is that of drawing the 
annuli on the celestial sphere, or on the corresponding 
projected plane. The algorithm we use based on LUTs 
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has proved to be more efficient than other methods we 
have studied, and this strategy is also employed by other 
groups [Hfij . 

A LUT is an array containing the list of borders of all 
the possible annuli, for a given value of v and /, clipped 
on the sky-patch we use. Therefore it contains the coordi- 
nates of the points belonging to the borders that intersect 
the sky-patch, in accordance to the tiling we use, together 
with information to take care of edge effects. As de- 
scribed in 1 2 ■''{!] , it turns out that the annuli are relatively 
insensitive to the value of the search frequency and, once 
a LUT has been constructed for a particular frequency, it 
can be re-used for a large number of neighboring frequen- 
cies thus allowing for computational savings. The value 
of v used to construct the LUTs corresponds to the aver- 
age velocity of the detector in the SSB frame during the 
30 minutes interval of the corresponding SFT. 

In fact, the code is further sped up by using partial 
Hough map derivatives (PHMDs) instead of the PHMs, 
in which only the borders of the annuli are indicated. A 
PHMD consists of only ones, zeros, and minus ones, in 
such a way that by integrating appropriately over the dif- 
ferent sky locations one recovers the corresponding PHM. 
This integration is performed at a later stage, and just 
once, after summing the appropriate PHMDs, to obtain 
the final Hough map. 

In the pipeline, we loop over frequency and spin-down 
values, taking care to update the set of PHMDs currently 
used, and checking the validity of the LUTs. As soon as 
the LUTs are no longer valid, the code recomputes them 
again together with the sky grid. Statistical analyses 
arc performed on the Hough maps in order to compress 
the output size. These include finding the maximum, 
minimum, mean and standard deviation of the number 
counts for each individual map, the parameters of the 
loudest event, and also of all the candidates above a cer- 
tain threshold. We also record the maximum number 
count per frequency analyzed, maximized over all spin- 
down values and sky locations, and a histogram of the 
number counts for each 1 Hz band. The schematic search 
pipeline is shown in Fig. [3J 

As a technical implementation detail, the search is 
performed by dividing the 200 Hz frequency band into 
smaller bands of 1 Hz and distributed using Condor |Hl| 
on a computer cluster. Each CPU analyzes a differ- 
ent 1 Hz band using the same pipeline (as described in 
Fig. |3J). The code itself takes care to read in the proper 
frequency band from the SFTs. This includes the search 
band plus an extra interval to accommodate for the max- 
imum Doppler shift, spin-down, and the block sized used 
by the running median. The analysis described here was 
carried out on the Merlin cluster at AEI j52j. The full- 
sky search for the entire S2 data from the three detectors 
distributed on 200 CPUs on Merlin lasted less than half 
a day. 

The software used in the analysis is available in 
the LIGO Scientific Collaboration's CVS archives at 
http : //www . lsc-group . phys . uwm . edu/cgi-bin/ cvs/ 




50 100 150 200 

Number count 




200 250 300 350 400 450 500 550 
Number count 




250 300 350 
Number count 



FIG. 4: Top: Graph of the LI number-count discrete proba- 
bility distribution: the solid curve corresponds to the number 
count distribution obtained for the band between 206-207 Hz, 
the dash-dot curve to the number count distribution obtained 
for the 343-344 Hz band, that contains violin modes, and in 
asterisks the theoretical expected binomial distribution for 
687 SFTs and a peak selection probability of 20%. Middle: 
the HI number-count distribution for 1761 SFTs. Bottom: 
the H2 number-count distribution for 1384 SFTs. 



viewevs . cgi/?cvsroot=lscsof t, together with a suite 
of test programs, especially for visualizing the Hough 
LUTs. The full search pipeline has also been validated 
by comparing the results with independently written 
code that implements a less efficient but conceptually 
simpler approach, i.e., for each point in parameter 
space (/o,/,n), it finds the corresponding pattern in 
the time-frequency plane and sums the corresponding 
selected frequency bins. 



D. Number counts from LI, HI and H2 

In the absence of a signal, the distribution of the Hough 
number count ideally is a binomial distribution. Environ- 
mental and instrumental noise sources can excite the op- 
tically sensed cavity length, or get into the output signal 
in some other way, and show up as spectral disturbances, 
such as lines. If no data conditioning is applied, line in- 
terference can produce an excess of number counts in the 
Hough maps and mask signals from a wide area in the 
sky. Fig. 0| shows the comparison of the theoretical bi- 
nomial distribution Eq. 1|12|) with the distributions that 
we obtain experimentally in two bands: 206-207 Hz and 
343-344 Hz. The first band contains very little spectral 
disturbances while the second band contains some violin 
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FIG. 5: Graph of the LI maximum number count 



400 



for 



every analyzed frequency fk, maximized over all spin-down 
values and sky locations. The top figure corresponds to the 
raw output from the Hough transform in which many out- 
liers are clearly visible. The bottom figure corresponds to the 
same data after vetoing the frequency bands contaminated by 
known instrumental noise. See App.|^|for details on outliers. 



modes. As shown in Fig. 21 the Hough number count 
follows the expected binomial distribution for the clean 
band while it diverges from the expected distribution in 
the presence of strong spectral disturbances, such as the 
violin modes in this case. We have verified good agree- 
ment in several different frequency bands that were free 
of strong spectral disturbances. 

The sources of the disturbances present in the S2 data 
are mostly understood. They consist of calibration lines, 
broad 60 Hz power line harmonics, multiples of 16 Hz due 
to the data acquisition system, and a number of mechan- 
ical resonances, as for example the violin modes of the 
mirror suspensions [36l |46| . The 60 Hz power lines are 
rather broad, with a width of about ±0.5 Hz, while the 
calibration lines and the 16 Hz data acquisition lines are 
confined to a single frequency bin. A frequency comb is 
also present in the data, having fundamental frequency 
at 36.867Hz for LI, 36.944Hz for HI and 36.975 Hz for 
H2, some of them accompanied with side lobes at about 
0.7 Hz, created by up-conversion of the pendulum modes 
of some core optics, but these were only present (or at 
least prominent) in HI and H2. The sources of these 
combs were synthesized oscillators used for phase modu- 
lation that were later replaced by crystal oscillators. In 
addition to the above disturbances, we also observe a 
large number of multiples of 0.25 Hz. While this comb of 
lines is strongly suspected to be instrumental, its physical 
origin is unknown. In Table [I] wc summarize the list of 
known spectral disturbances in the three interferometers 
during the S2 run. 



LI 


HI 


H2 


/ (Hz) 


A/ (Hz) 


/ (Hz) 


A/ (Hz) 


/ (Hz) 


A/ (Hz) 


0.250* 


0.00 


0.250* 


0.00 


0.250* 


0.00 


16.000* 


0.00 


16.000* 


0.00 


16.000* 


0.00 


60.000* 


1.00 


60.000* 


1.00 


60.000* 


1.00 


221.200 


0.01 


221.665 


0.01 


221.850 


0.02 


258.080 


0.04 


257.875 


0.02 


258.830 


0.01 


294.935 


0.01 


258.610 


0.02 


295.070 


0.02 


331.810 


0.01 


259.340 


0.02 


295.670 


0.02 


345.000 


6.00 


294.820 


0.01 


295.800 


0.04 


368.670 


0.02 


295.560 


0.00 


295.930 


0.02 






296.300 


0.00 


296.530 


0.02 






331.790 


0.00 


323.300 


0.00 






332.490 


0.00 


323.870 


0.04 






333.200 


0.00 


324.000 


0.04 






335.780 


0.14 


324.130 


0.04 






336.062 


0.00 


324.700 


0.00 






339.000 


0.02 


332.800 


0.00 






339.720 


0.02 


335.120 


0.02 






345.000 


6.00 


335.590 


0.02 






365.500 


0.02 


341.615 


0.01 






368.690 


0.00 


346.500 


9.00 






369.430 


0.00 


349.202 


0.00 






370.170 


0.01 







★ indicates all higher harmonics are present 

TABLE I: List of known spectral disturbances in the three 
interferometers during the S2 run used as a frequency veto 
in the 200-400 Hz band. / refers to the central frequency 
and A/ to the full width of the lines. Lines denoted with 
A/ = 0.0 Hz are those in which the line width is much smaller 
than the associated maximum Doppler broadening of the line. 
This ranges from ~ 0.04 Hz at a frequency of 200 Hz up to 
~ 0.08 Hz at 400 Hz. 
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FIG. 6: Graph of the HI maximum number count n* k versus 
frequency fk as in Fig. 
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FIG. 7: Graph of the H2 maximum number count nl versus 
frequency fk as in Fig. [K] 
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FIG. 8: Histograms of the maximum Hough number count 
n k for the three detectors. The light curve corresponds to the 
raw output and the dark thick curve corresponds to the same 
output after vetoing the contaminated frequencies. 



After analyzing all the data we discard the number 
counts from all those frequencies that could be con- 
taminated by a known instrumental artifact. Thus, we 
exclude every frequency bin which is affected by the 
spectral disturbances including the maximum possible 
Dopplcr broadening of these lines; thus, for a known 









Before Veto 


After Veto 




«.th(o20o) 


J1th( Q 40u) 


(nl) std«) 


(nl) std(n£) 


LI 


188.8 


191.7 


194.8 32.3 


191.4 10.2 


HI 


435.8 


440.3 


452.0 94.0 


439.8 29.5 


H2 


350.6 


354.6 


360.6 72.3 


353.2 11.3 



TABLE II: Comparison between the statistics of the max- 
imum Hough number count nl (before and after the fre- 
quency veto) and n th (a) at different false alarm rates 0:200 = 
9.5 x 10 -7 and 0:400 = 2.4 x 10~ 7 for the three detectors. 



spectral disturbance at a frequency / and width A/, we 
exclude a frequency range ±(vf/c + A//2) around the 
frequency /. We also exclude from our analysis the fre- 
quency band 342-348 Hz for LI and HI and 342-351 Hz 
for H2 because they contain many violin modes. The net 
effect of this vetoing strategy is that we consider only 
67.1% for LI, 66.8% for HI and 65.6% for H2 of the full 
200 Hz range. 

Figures [SJ and [7] show the maximum Hough num- 
ber count n* k obtained in each of the 360000 different 
frequency bins fk analyzed, maximized over all sky lo- 
cations and spin-down values for the three detectors. In 
those figures we compare the results of the search be- 
fore and after applying the frequency veto, showing how 
the spectral disturbances impact on the number counts. 
This can be seen more clearly in Fig. [S] where we plot 
the histograms of these maximum Hough number counts 
n* k before and after applying the frequency veto. 

These values of n* k obtained after removing the known 
outliers can be easily seen to be consistent with what 
we would expect for Gaussian stationary noise. As- 
suming that the maximum occurs only once, the ex- 
pected value of rif. should be consistent with a false 
alarm rate of l/nz bin (/fc), where m bin (/fe) is the total 
number of templates at a frequency fk- rn bin (f k ) is fre- 
quency dependent, ranging from ~ 1.1 x 10 6 at 200 Hz 
to — 4.2 x 10 6 at 400 Hz. Thus, n£ should be consistent 
with a false alarm rate of 0:400 = 2.4 x 10~ 7 at 400 Hz, 
up to 0:200 = 9.5 x 10~ 7 at 200 Hz. In the case of Gaus- 
sian stationary noise, the expectation value of n k should 
therefore be similar to n th (defined by Eq. (|16fl ) for a 
false alarm a = l/m bin (fk). Table ITT1 compares the mean 
(n k ) of the maximum Hough number count before and 
after vetoing, with n th (a) at different false alarm rates. 
After vetoing, we observe that (n* k ) lies within the in- 
terval (n th (0:200); ^th(a40o))j and the standard deviation 
std(n£) is also greatly reduced. This indicates the con- 
sistency of the observed values of n k with ideal noise. 

As can be seen in Fig.[3|Sland|3 a few outliers remain 
after applying the frequency veto described above. The 
fact that these outliers have such small false alarm prob- 
abilities makes it very unlikely they were drawn from a 
parent Gaussian distribution. However they could also 
be due to spectral disturbances (line noise) that mimic 
the time-frequency evolution of a pulsar for a certain lo- 
cation in the sky. If these outliers are due to gravita- 
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tional signals they should show up in the different de- 
tectors. By performing a simple coincidence analysis in 
frequency, the only candidate that remains is the one at 
210.360 Hz. But this has been ruled out since it seems 
to be associated to multiples of 70.120 Hz produced by 
a VME (VERSA module Eurocard) controller hardware 
used during S2. Since the other outliers are not coinci- 
dent among the three detectors, there is no evidence for a 
detection. We refer the interested reader to appendix 151 
for further details. 

As explained in Sec. [fl the ultimate goal for wide pa- 
rameter space searches for continuous signals over large 
data sets is to employ hierarchical schemes which alter- 
nate coherent and semi-coherent techniques. The Hough 
search would then be used to select candidates in the 
parameter space to be followed up. The way in which 
those candidates would be selected is the following: fix 
the number of candidates to follow up, determine the 
false alarm rate, and set the corresponding threshold on 
the Hough number count. Not all the candidates selected 
in this manner will correspond to real gravitational wave 
signals, but they will point to interesting areas in param- 
eter space. 

The analysis presented here is a very important step 
forward in this direction. However, given the limited sen- 
sitivity its relevance mostly rests in the demonstration of 
this analysis technique on real data. In what follows we 
will thus concentrate on setting upper limits on the am- 
plitude ho of continuous gravitational waves emitted at 
different frequencies. 



VII. UPPER LIMITS 

We use a frcqucntist method to set upper limits on the 
amplitude ho of the gravitational wave signal. Our upper 
limits refer to a hypothetical population of isolated spin- 
ning neutron stars which are uniformly distributed in the 
sky and have a spin-down rate / uniformly distributed in 
the range (— f max ,0). We also assume uniform distribu- 
tions for the parameters cost £ [—1,1], ip £ [0, 2w], and 
$o G [0, 27r]. As before, the frequency range considered 
is 200-400 Hz. 

The upper limits on ho emitted at different frequen- 
cies are based on the highest number count, the loudest 
event, registered over the entire sky and spin-down range 
at that frequency. Furthermore, we choose to set upper 
limits not on each single frequency but on a set of fre- 
quency values lying within the same 1 Hz band and thus 
are based on the highest number count in each frequency 
band. In every 1 Hz band the loudest event is selected 
excluding all the vetoed frequencies of Table [I] 

Let n* be the loudest number count measured from 
the data. The upper limit h§ on the gravitational wave 
amplitude, corresponding to a confidence level C, is the 
value such that had there existed in the data a real signal 
with an amplitude greater than or equal to /iq% then in 
an ensemble of identical experiments with different re- 
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FIG. 9: Histograms of the Hough number count distribution 
(in arbitrary units) for LI using 1000 injected signals ran- 
domly distributed over the whole sky within the band 200- 
201 Hz for different ho values. The largest number count for 
the search in that band was 202. The confidence level asso- 
ciated with the different ho values are: 0.1% for 1.0 x 10~ 23 , 
30.5% for 2.0 x 10 -23 , 87.0% for 4.0 x 10~ 23 , and 1 for 
1.0 x 10~ 22 . 



alizations of the noise, some fraction C of trials would 
yield a maximum number count equal to or bigger than 
n*. The upper limit hg corresponding to a confidence 
level C is thus defined as the solution to this equation: 



n 



Prob(n > n*\hg) = £ p(n\hg) = C, 



(23) 



hl 5% denotes the 95% 



where p(n\ho) is the number count distribution in the 
presence of a signal with amplitude ho and averaged over 
all the other parameters; note that p(n\ho) is different 
from the distribution p(n\h) discussed before Eq. I|12|) 
which was relevant for a targeted search when all signal 
parameters are known. We choose to set upper limits at 
a confidence level of C = 95% 
confidence upper limit value. 

Given the value of n* , if the distribution p(n\ho) were 
known, it would be a simple matter to solve Eq. I|23|) 
for h§ ■ In the absence of any signal, this distribution 
is indeed just a binomial, and as exemplified in Fig. 0] 
this is also what is observed in practice. If a signal were 
present, the distribution may not be sufficiently close to 
binomial because the quantity A defined in Eq. l|14(l varies 
across the SFTs due to non-stationarity in the noise and 
the amplitude modulation of the signal for different sky 
locations and pulsar orientations. In addition, now we 
must also consider the random mismatch between the 
signal and template (in relation to the parameter space 
resolution used in the analysis) which causes an addi- 
tional reduction in the effective SNR for the template. 
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FIG. 10: Pipeline description of the Monte-Carlo simulation 
to determine the upper limit values hg 5% . We inject randomly 
generated fake pulsar signals with fixed amplitude (the other 
parameters are drawn from a suitable uniform distribution) 
into the real data and measure the value of ho required to 
reach a 95% confidence level. 



For these reasons, we measure p(n\ho) by means of a sc- 
ries of software injections of fake signals in the real data. 
Fig. |5| shows four distributions for different ho values ob- 
tained by Monte Carlo simulations. While for low signal 
amplitudes the distribution is close to the ideal binomial 
one, the distribution diverges from it at higher ampli- 
tudes, thereby illustrating the complexity of the number 
count distribution for sufficient large ho- 

Our strategy for calculating the 95% upper limits is to 
find p(n\ho) for a wide range of ho values, then to get the 
corresponding confidence levels C(ho), and find the two 
values of ho which enclose the 95% confidence level. The 
95% upper limit is approximated by a linear interpolation 
between these values. We then refine this reduced range 
of ho values by further Monte-Carlo simulations until the 
desired accuracy is reached. This is described in Fig. El 

The parameters of the fake injected signals arc drawn 
from the population described above, and we ensure that 
the frequency does not lie within the excluded bands. 
The data with the injections are searched using the search 
pipeline used for the actual analysis. For computational 
efficiency, for each injected signal we find the number 
count using only the 16 nearest templates, and choose 
the one yielding the largest number count. The SFT data 
in the different frequency bins in a 1 Hz band (of order 
TV x 1800 bins that get combined differently for the differ- 
ent time-frequency patterns) can be considered as differ- 
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FIG. 11: Confidence level as a function of the signal ampli- 
tude ho for different Monte Carlo simulations for LI within 
the band 200-201 Hz. The solid line corresponds to a sim- 
ulation using 100000 injected signals. The two dashed lines 
correspond to 5000 injected signals. The dash-dotted line cor- 
responds to 3000 signal injections. The /ig 5% upper limit for 
this band using and these two simulations with 5000 injections 
has a maximum absolute error of the order of 0.02 x 10~ 2,i 
corresponding to a relative error smaller than 0.5%. In the 
case of using only 3000 injections the error increases to the 
1.5% level. 



ent realizations of the same random process. Therefore 
it is reasonable to assume that the normalized histogram 
of the largest number count represents the discrete prob- 
ability distribution p(n\ho). 

The most computationally intensive part of this 
Monte-Carlo scheme is the generation of the artificial sig- 
nals. The computational costs can be greatly reduced by 
estimating p{n\h<y} for different ho values in one go: for 
each individual artificial signal, we generate a set of SFTs 
containing only one noiseless signal of a given amplitude. 
These SFTs can be scaled by an appropriate numerical 
constant to obtain a set of SFTs containing signals with 
different amplitudes, which are then added to the noise 
SFTs. The disadvantage of doing this is that the differ- 
ent signals obtained by rescaling the amplitude this way 
are not statistically independent since all the other signal 
parameters are identical. We must ensure that we have 
a sufficiently large number of trials so that the error in 
the final upper limit is sufficiently small. 

We have found empirically that 5000 injections per 
band are sufficient to get upper limits accurate to within 
3%; sec Fig. El d and US Fig. HP shows the confidence 
level as a function of the signal amplitude ho for LI data 
within the band 200-201 Hz. The solid line corresponds 
to our most accurate simulations using 100,000 injections. 
The two dashed lines correspond to two different simula- 
tions both using 5000 injected signals and the dotted line 
corresponds to simulation with 3000 injections. In each 
case, confidence levels for different ho values arc calcu- 
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FIG. 12: Confidence level as a function of the signal ampli- 
tude ho for different Monte Carlo simulations for LI within 
the band 200-201 Hz. The solid thick line corresponds to a 
simulation using 100000 injected signals. The other 20 lines 
correspond to simulations with 5000 injected signals. The 
hg 5% upper limit for this band using 5000 injections has a 
maximum absolute error of the order of 0.1 x 10 -23 corre- 
sponding to a relative error of 2.2%. 




FIG. 13: Confidence level as a function of the signal ampli- 
tude ho for different Monte Carlo simulations for HI within 
the band 259-260 Hz. The solid thick line corresponds to a 
simulation using 100000 injected signals. The other 20 lines 
correspond to simulations with 5000 injected signals. The 
ho 5% upper limit for this band using 5000 injections has a 
maximum absolute error of the order of 0.1 x 10 -23 corre- 
sponding to a relative error of 2.0%. 



lated by simply rescaling the signal as described above. 
This means that all the points in each individual curve 
in Fig. ^2 are statistically biased in a similar way, and 
this explains why the curves in this figure do not inter- 
sect. To estimate the error in the 95% upper limit, we 
generate several of these curves for a fixed number of in- 
jections and we measure the error in /jq 5 ^ with respect 
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FIG. 14: Measured Hough number count probability distri- 
bution for p(n|/iQ njcot ) from the Monte Carlo simulations using 
100,000 injected signals, (a) LI within the band 200-201 Hz, 
with a h n ' e<:t of 4.422 x 10~ 23 corresponding to a confidence 
level of 94.95%. The n* value for this band was 202. (b) HI 
within the band 259-260 Hz, with a /ij," icct of 4.883 x 10" 23 
corresponding to a confidence level of 95.04%. The n* value 
for this band was 455. (c) H2 within the band 258-259 Hz, 
with a /iQ njcct of 8.328 x 10 -23 corresponding to a confidence 
level of 95.02%. The n* value for this band was 367. 



to the accurate reference value obtained using 100,000 
injections. For this particular band and detector, we find 
that for 5000 injections, the error in the upper limit is 
at most 0.1 x 10 -23 , corresponding to a relative error of 
2.2%. This experiment has been repeated for several 1 Hz 
bands and many simulations, and in all of them 5000 
injections were enough to ensure an accuracy of better 
than 3% in the /i^ 5% . We also have found that with 5000 
injections, using amplitudes equal to the most accurate 
upper limit h^ 5% obtained from 100,000 injections, the 
confidence level are within 94.5 — 95.5%. See Fig. 1121 and 

ED 

Fig. [21 shows the probability distribution for 
p(n|ft.Q joct ) for LI, HI and H2 in a 1 Hz band measured 
from 100,000 randomly injected signals, with a signal am- 
plitude /iQ Jcct close to the 95% confidence level; note that 
these /io JBCt values are not the 95% upper limits because 
they correspond to different confidence levels. This illus- 
trates, yet again, that the true number count distribu- 
tions are far from the ideal binomial distribution. It is 
also interesting to notice that, among all the templates 
used to analyze each injection, the nearest template, in 
the normal Euclidean sense, corresponds to the template 
providing the highest number count only 18% of the time. 
This effect is due to the noise contribution and the fact 
that the matching-amplitudc, described by the paramc- 
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FIG. 15: The 95% confidence upper limits on ho over the 
whole sky and different spin-down values in 1 Hz bands. 



ter space metric is highly anisotropic |53j . 

The 95% confidence upper limits on ho for each 1 Hz 
band using all the data from the S2 run are shown in 
Fig. As expected, the results are very similar for the 
LI and HI interferometers, but significantly worse for H2. 
The typical upper limits in this frequency range for LI 
and HI are mostly between ~ 4-9 x 10~ 23 , typically bet- 
ter at lower frequencies. The most stringent upper limit 
for LI is 4.43 x 1CT 23 which occurs within 200-201 Hz, 
largely reflecting the lower noise floor around 200 Hz. For 
HI it is 4.88 x 10 - 23 , which occurs in the frequency range 
259-260 Hz, and for H2 it is 8.32 x 10" 23 , which occurs in 
the frequency range 258-259 Hz. The values of the most 
stringent upper limits on /ig 5% have been obtained using 
100,000 injections in the most sensitive 1 Hz bands. The 
upper limits are significantly worse in bands lying near 
the known spectral disturbances, especially near the vi- 
olin modes. 

In Table II I II we summarize the best upper limits on 
hf % and we compare them with the theoretical values 
h° Q xp we would expect for a directed search using a per- 
fectly matched template, as given by Eq. l(TT}) . Here we 
take a false dismissal rate of 5% and the false alarm rate 
associated to the loudest number count in that band. In 



those three bands the ratio /ip 5% //ig xp is about 1.8. 

These hf % results are also about a factor of 2.6 worse 
than those predicted by Fig. ^ which corresponds to a 
directed search using a perfectly matched template, and 
with a false alarm rate of 1% and a false dismissal rate 
of 10%. Of these 2.6, a factor ~ 1.5 is due to the use 
of different values of the false alarm and false dismissal 
rate, and a factor ~ 1.8 because the p(n\ho) distribution 
does not correspond to the ideal binomial one for values 
of ho distinct from zero. 

From the upper limits on ho, using Eq. (JSJ, we can 
derive the distance covered by this search. This is shown 
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FIG. 16: Astrophysical reach covered by the search. The 
curves show the estimated distance d out to which signals 
from isolated gravitational wave pulsars could be detected in 
our S2 data set derived from the upper limits on ho. This 
plot assumes e = 10 -6 and I zz = 10 45 gcm 2 . 



in Fig. 1161 assuming e = 1Q~ 6 and I zz = 10 45 gcm 2 . The 
maximum reach is 2.60 pc for HI which occurs within 
395-396 Hz. For LI the maximum reach is 2.15 pc and 
1.62 pc for H2. This value of e = 10~ 6 corresponds to the 
maximum expected ellipticity for a regular neutron star, 
but ellipticities from more exotic alternatives to neutron 
stars may be larger |10|. e.g., solid strange quark stars 
for which e max 10 -4 . Therefore the astrophysical reach 
for these exotic stars could be better by a factor 100. 



VIII. HARDWARE INJECTIONS 



Two artificial pulsar signals, based on the waveforms 
given in Eqs. 0, and (jfJJ) were injected into 

all three LIGO interferometers for a period of 12 hours 
towards the end of the S2 run. These injections were 



designed to give an end-to-end validation of the search 
pipeline. The waveforms were added digitally into the 
interferometer Length Sensing and Control System (re- 
sponsible for maintaining a given interferometer on res- 
onance), resulting in a differential length dither in the 
optical cavities of the detector. We denote the two pul- 
sars PI and P2; their parameters are given in Table llVl 
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Detector 


Best hf % 


/band (HZ) 


<S„(/ b „ d )> (Hz- 1 ) 


N 


n 


a" 




miHz 


5 oxp 


Kr 




LI 


4.43 x 1(T 23 


200-201 


1.77 x 10" 43 


687 


202 


8.94 x 1(T 


10 


1.88 x 10 9 


5.4171 


2.41 x 10" 


23 


HI 


4.88 x 1(T 23 


259-260 


3.53 x 10~ 43 


1761 


455 


1.77 x 10" 


-9 


3.11 x 10° 


5.3381 


2.67 x 10" 


_>:! 


H2 


8.32 x 10" 23 


258-259 


1.01 x 10~ 42 


1384 


367 


2.25 x 10" 


-9 


3.11 x 10 9 


5.3098 


4.78 x 10" 


23 



TABLE III: Best all-sky upper limits obtained on the strength of gravitational waves from isolated neutron stars. The 
values have been obtained using 100,000 injections in the best 1 Hz bands /band for the three detectors. {Sn(f bB , ad )) is the 
average value of noise in that 1 Hz band excluding the vetoed frequencies, N the number of SFTs available for the entire 
S2 run, n* is the loudest number-count measured from the data in that band, a* is the corresponding false alarm assuming 
Gaussian stationary noise derived from Eq. 11611 . mi Hz the number of templates analyzed in that 1 Hz band, the quantify 5 exp 
is defined by Eq. 11811 using the values a = a* and j3 = 0.05. h^ p is the theoretical expected upper limit from searches with 
perfectly matched templates assuming the ideal binomial distribution for p(n\ho) defined by Eq. (U), therefore ignoring also 
the effects of the different sensitivity at different sky locations and pulsar orientations. 



fa (Hz) 
/ (Hz-s- 1 ) 
RA (rad) 
Dec (rad) 
tj) (rad) 
cos i 
$0 (rad) 
To (sec) 
ho 



PI 

1279.123456789012 


5.1471621319 
0.3766960246 







733967667.026112310 
2 x 10 -21 



P2 

1288.901234567890123 

-10 -8 
2.345678901234567890 
1.23456789012345 




733967751.522490380 
2 x 10~ 21 



TABLE IV: Parameters of the two hardware injected pulsars. 
See Eqs. {J, @, |@}, and @ for the definition of the param- 
eters. RA and Dec are the right ascension and declination in 
equatorial coordinates. To is the GPS time in the SSB frame 
in which the signal parameters are defined, ho is the ampli- 
tude of the signal according to the strain calibration used at 
the time of the injections. 



The data corresponding to the injection period have 
been analyzed using the Hough transform and the same 
search code as described in Sec. IVI CI As before, the 
input data consists of 30 minutes long SFTs. The num- 
ber of SFTs available are 14 for LI, 17 for HI and 13 
As in 



for H2. As in Sec. IVI Bl the frequency resolution is 
1 /1800 s, the sky resolution is given by equation (|2*2*)l , and 

Sf = 1 / (T ob8 T coh ) . Since the total effective observation 
time is somewhat different for the three detectors, we get 

1.77024 xl0" 8 Hz- 
1 for H2. As be- 



Sf 



-2.28624x 10~ 8 Hz-s" 1 for LI, 



s~ 4 for HI, and -1.93533 x 10~ 8 Hz-s 
fore, for each intrinsic frequency we analyze 11 different 
spin-down values. The portion of sky analyzed was of 
0.5 radians x 0.5 radians around the location of the two 
injected signals. 

Fig. El shows the Hough maps corresponding to the 
nearest frequency and spin-down values to the injected 
ones. Although the presence of the signal is clearly visi- 
ble, it is apparent that 12 hours of integration time is 
not enough to identify the location of the source us- 
ing the Hough transform. In particular if one looks at 
the Hough maps at mismatched frequencies and spin- 
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FIG. 17: Hough maps for the hardware injected signals. The 
left maps correspond to PI and the right maps to P2. From 
top to bottom, the left maps correspond to LI, HI and H2, for 
a frequency of 1279.123333 Hz and zero spin-down; the right 
maps to a frequency of 1288.901111 Hz and zero, —1.77024 x 
10- 8 Hz-s _1 , and -1.93533 x 10 -8 Hz-s -1 spin-down values 
respectively. The location of the injected signals are close to 
the centers of these subplots. 



down values, one can still identify annuli with very high 
number-counts, but appearing with a mismatched sky 
location. 

For P2, the Hough maps corresponding to the closest 
values of frequency and spin-down contain the maximum 
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Detector / band (Hz) hf 



LI 


200-201 


4.43 x 10" 


- 23 


HI 


259-260 


4.88 x 10" 


-23 


H2 


258-259 


8.32 x 10" 


- 23 



TABLE V: Best all-sky upper limits obtained on the strength 
of gravitational waves from isolated neutron stars. The hg 5% 
values have been obtained using 100,000 injections in the best 
1 Hz bands. 



number count at the correct sky location. These maxi- 
mum number counts are 12 for LI, 17 for HI, and 11 for 
H2. Notice that for LI and H2 these numbers are smaller 
than the number of SFTs used. 

For PI, the closest template to the signal parameters 
has a number count of 13 for LI, 12 for HI and 8 for 
H2. The maximum number counts obtained in the search 
were 13 for LI, but 14 for HI, and 10 for H2. Those 
higher number counts occurred for several templates with 
a larger mismatch, for example at 1279.123333 Hz and 
-5.31073 x lO^Hz-s- 1 for HI, and 1279.132222Hz and 
-1.35473 x lO^Hz-s -1 for H2. This is not surprising 
because we only compute the Hough maps at the Fourier 
frequencies n x l/T coh . In any case, both pulsar sig- 
nals are unambiguously detected because these number 
counts are much bigger then the expected average num- 
ber counts for pure noise. 



IX. CONCLUSIONS 

In this paper we use the Hough transform to search 
for periodic gravitational- wave signals. This is a semi- 
coherent sub-optimal method. Its virtues are computa- 
tional efficiency and robustness. The search pipeline was 
validated using a series of unit tests and comparisons 
with independently written code. We then applied this 
method to analyze data from the second science run of 
all three LIGO interferometers. We also validated the 
search pipeline by analyzing data from times when two 
artificial pulsar signals were physically injected into the 
interferometer hardware itself. We show in this paper 
that the injected signals were clearly detected. 

Our final results are all-sky upper limits on the gravita- 
tional wave amplitude for a set of frequency bands. The 
best upper limits that we obtained for the three inter- 
ferometers are given in Table ]V\ The overall best upper 
limit is 4.43 x 10 -23 . We searched the 200 - 400 Hz band 
and the spin-down space / max < 1.1 x 10 _9 Hz-s _1 , and 
no vetoes were applied except for the list of ignored fre- 
quency bands that contain instrument line artifacts. Our 
best upper limit is 26 times worse than the best upper 
limit obtained for a targeted coherent search using the 
same data. This was an upper limit of 1.7 x 10 -24 [l6[ . 
achieved for PSR J1910-5959D. This is to be expected 
because we have performed not a targeted, but a wide 
parameter space search. If we were to use the optimal 



J-"-statistic method to perform a hypothetical search over 
the same parameter space region as the Hough search in 
this paper, the number of templates required would be 
much larger for the same observation time: ~ 10 19 [2(j 
instead of ~ 10 12 . Thus, we would have to set a lower 
false alarm rate for this hypothetical search, and in the 
end, the sensitivity turns out to be roughly comparable 
to that of the Hough search. Note also that this hy- 
pothetical search is not computationally feasible for the 
foreseeable future. 

We can search for smaller signals either by increasing 
the number of coherent segments, TV, or by increasing 
the coherent time baseline, T coh . Since the number of 
segments is determined by the length of the data set, for 
a given amount of data one wants to make T coh as large 
as possible. However, for the search pipeline presented 
in this paper, increasing T coh was not possible due to the 
restriction on its value mentioned in Sec. [V] This will 
be overcome by demodulating each short segment, taking 
into account the frequency and amplitude modulation of 
the signal. The optimal method will be to calculate the 
^-"-statistic ^3 f° r each segment. The time- frequency 
pattern will then no longer be given by (JSJ) and i|10|) but 
by the master equation given in reference |23l |. 

Since the wide parameter-space search for periodic 
gravitational-wave signals is computationally limited, 
there is also a limit on the maximum T coh that can be 
used, given finite computing resources. Thus, a hierar- 
chical strategy that combines fully-coherent and semi- 
coherent methods will be needed to achieve optimal re- 
sults 0, ■ Our goal is to use the Hough transform as 
part of such a strategy. This is work in progress and the 
results will be presented elsewhere. 
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FIG. 18: Value of the running median bias C as a function 
of the window size w, C approaches log(2.0) for large w. 



APPENDIX A: 



THE BIAS IN THE RUNNING 
MEDIAN 



We are using a running median to estimate the noise 
floor in our SFTs. Thus, the value of S n at a particular 
frequency bin can be estimated from the median of \xk\ 2 
in w frequency bins around the frequency bin, where w is 
an integer and represents the window-size of the running 
median. The reason for using a running median is to 
minimize the effect of large spectral disturbances, and to 
ensure that the presence of any pulsar signals will not 
bias our estimation of S n . To carry this out in practice, 
we would like to know how the median can be used as an 
estimator of the mean. In this appendix, we answer this 
question assuming that the noise is Gaussian, so that the 
power is distributed exponentially 

Let £ be a random variable with probability distribu- 
tion f(x). Let F(x) denote the cumulative distribution 
function: 



F(x) 



f(x')dx' . 



(Al) 



Let us draw w samples from this distribution and arrange 
them in increasing order: x n (n = l...w). Define an 
integer k which is (w — l)/2 when w is odd, and w/2 
when w is even. We define the median as 



*(«)) _ J Xk+\ when w is odd , 

I \{xk + Xk+i) when w is even. 



(A2) 



Consider first the case when w is odd. The distribution 
of can be found as follows: lies within the range 
(x, x + dx) when k values are less than x, w — k — l values 
are greater than x + dx, and one value is in the range 



(x, x + dx). The probability density for is thus 



g w {x)=\ W k (w - k)[F(x)] k [l - F(x)] 1 



-fc-i 



fix). 



(A3) 

When w is very large, it can be shown that the distribu- 
tion g w (x) approaches a Gaussian whose mean is equal 
to the population median £, and whose variance is pro- 
portional to 1/ y/w [54| ■ 

For Gaussian noise, the normalized power p follows 
the exponential distribution, i.e. f(x) = e~ x for x > 
and f(x) = for x < 0. The mean is unity, therefore 
is an unbiased estimator of the mean, where 
E(C (tu) ) is the expectation value of C {w) - 



E(C W ) = / 

J —1 



xg w (x)dx . 



(A4) 



We can explicitly calculate E^™') for this case and the 
answer turns out to be given by a truncated alternating 
harmonic series: 



e(c m ) = E 



(_l)i+i 



J 



(A5) 



For very large w, ~E(£( W " > ) approaches ln(2), which is pre- 
cisely the population median. For w = 1 it is just unity, 
which makes sense because in this case the median is 
equal to the mean. For finite w, E^™)) is somewhat 
larger than ln(2) and for a window size of w = 101, which 
is what is used in the actual search, it is 0.698073. This 
is to be compared with ln(2) = 0.693147, a difference of 
about 0.7%. 

When w is even, the distribution of C^ w ' is given (up 
to a factor of 2) by the convolution of the distributions 
of Xk and Xk+i- However, we are interested only in the 
expectation value 



E(C^) = \[E(x k ) 



E(z fc+ i 



(A6) 



The expectation value of Xk+i is calculated as above, us- 
ing the distribution (|A3(I . while the distribution of x k 
is obtained by replacing k with k — 1 in i|A3() . It turns 
out that E(C (u,) ) = E(C ( "'~ :L) ) when w is even. Thus 
E(C (2) ) = E(C (1) ), E(C (4) ) = E(C (3) ) and so on. Fig- 
urelUplots E(C (u,) ) for all values of w from 1 to 200. 



APPENDIX B: THE NUMBER COUNT 
OUTLIERS 

This appendix contains a discussion about the outliers 
in the Hough number counts that are strongly suspected 
to be instrumental artifacts but that we were not able to 
definitely identify as such. 

After applying the frequency veto described in Sec. 
II Dl we focus our attention on those candidates with 
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Detector Frequency (Hz) 


Number count 


LI 


210.36 


268 




234.50 


224 




281.35 


218 




329.34 


626 




335.62 


219 


HI 


210.36 


596 




212.26 


519 




244.14 


507 




249.70 


746 




280.48 


949 




329.58 


1510 




329.78 


1227 




348.45 


482 




350.60 


1423 


H2 


202.18 


402 




203.23 


395 




210.36 


443 




298.81 


394 




329.69 


867 




387.05 


400 




389.40 


391 



TABLE VI: List of outliers present in the Hough maps after 
applying the known instrumental frequency veto, for a false 
alarm rate of 10 -13 . This corresponds to a threshold in the 
number count of 216 for LI, 480 for HI and 390 for H2. For 
each outlier we quote the central frequency and the maximum 
number count. The triple coincidence at 210.36 Hz is a har- 
monic of the 70.12 Hz spectral disturbance described in the 
text. 



a false alarm rate a (for a single detector) smaller than 
10~ 13 . This corresponds to a threshold on the number 
count of 216 for LI, 480 for HI and 390 for H2. Since 
the total number of templates analyzed in this 200 Hz 
search band is roughly 10 12 , the probability of getting 
one candidate above that threshold over the full search is 
approximately 10% in each detector. All the candidates 
that satisfy such condition tend to cluster around a few 
frequencies that arc listed in Table IVII These are the 
so-called outliers that were present in Fig. [3] and [7| 

If these outliers are due to gravitational signals they 
should show up in the different detectors. By performing 
a simple coincidence analysis in frequency, the only can- 
didate that remains, at this false alarm level, is the one at 
210.36 Hz. The reader should notice that 210.36 Hz cor- 
responds to 3 x 70.12 Hz. In the HI data, there are also 
excess of number counts at 280.480 Hz and 350.600 Hz, 
corresponding to 4 x 70.12 and 5 x 70.12 respectively 
These 70.120 Hz multiples together with the 244.14 Hz 
line were detected in association with a VME (VERSA 
module Eurocard) controller hardware used during S2. 
However, since the data acquisition system architecture 
has changed since S2 the coupling mechanism cannot be 
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FIG. 19: The square root of the average value of S„ using the 
entire S2 LI data set analyzed. The four graphs correspond to 
the frequencies where outliers were present in the Hough maps 
after applying the known instrumental veto. They correspond 
to 210.36 Hz, 234.50 Hz, 281.35 Hz, 329.34 Hz and 335.62 Hz. 



proven. 

Figs. 1191 1201 and !21l show how the the outliers listed in 
Tablc lVll stand well above the background noise spectrum 
level S n , when this is estimated from the entire run. We 
believe they all arise from instrumental or environmental 
artifacts. However we are not able to determine in a con- 
clusive manner their physical cause. For this reason, the 
potential contaminated frequency bands have not been 
vetoed when setting the upper limits. 
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FIG. 20: The square root of the average value of S„ 
using the entire S2 HI data set analyzed. The graphs 
correspond to zooms in the frequencies where outliers 
were present in the Hough maps. They correspond 
to 210.36 Hz, 212.26 Hz, 244.14 Hz, 249.70 Hz, 280.48 Hz, 
329.58 Hz, 329.78 Hz, 348.45 Hz and 350.60 Hz. 
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